knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(beepr)


source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = TRUE
report_ordinal = FALSE

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() beepr::beep(sound = 5)
)
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

# Importantly, we do not recode pushing here. 
df_double <- prepare_data(df, recode_pushing = FALSE, use_mi = use_mi)[[1]]

Constructing scales reshaping data (4field) centering data within and between

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'ma[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'Daily received persuasion target -> target', 
    'Daily received persuasion target -> agent', 
    'Daily received pressure target -> target', 
    'Daily received pressure target -> agent', 
    'Daily received pushing target -> target', 
    'Daily received pushing target -> agent', 
    'Day', 
    'Daily weartime',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'Mean received persuasion target -> target',
    'Mean received persuasion target -> agent',
    'Mean received pressure target -> target',
    'Mean received pressure target -> agent',
    'Mean received pushing target -> target',
    'Mean received pushing target -> agent',
    'Mean weartime'
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Daily received persuasion target -> target)', 
  'sd(Daily received persuasion target -> agent)', 
  'sd(Daily received pressure target -> target)', 
  'sd(Daily received pressure target -> agent)', 
  'sd(Daily received pushing target -> target)', 
  'sd(Daily received pushing target -> agent)', 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]',
  'ma[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'ma[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,9),
  "Between-Person Effects" = c(10,16),
  "Random Effects" = c(17, 23), 
  "Additional Parameters" = c(24,29)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,9+5),
  "Between-Person Effects" = c(10+5,16+5),
  "Random Effects" = c(17+5, 23+5), 
  "Additional Parameters" = c(24+5,29+6)
  )
HURDLE MODELS
# For indistinguishable Dyads
model_rows_fixed_hu <- c(
    'Intercept', 
    'hu_Intercept',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    
    # HURDLE MODEL
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'hu_persuasion_self_cw', 
    'hu_persuasion_partner_cw', 
    'hu_pressure_self_cw', 
    'hu_pressure_partner_cw', 
    'hu_pushing_self_cw', 
    'hu_pushing_partner_cw', 
    'hu_day', 
    'hu_weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'hu_persuasion_self_cb',
    'hu_persuasion_partner_cb',
    'hu_pressure_self_cb',
    'hu_pressure_partner_cb',
    'hu_pushing_self_cb',
    'hu_pushing_partner_cb',
    'hu_weartime_self_cb'
  )

model_rows_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(hu_Intercept)',
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # HURDLE
  'sd(hu_persuasion_self_cw)',
  'sd(hu_persuasion_partner_cw)',
  'sd(hu_pressure_self_cw)',
  'sd(hu_pressure_partner_cw)',
  'sd(hu_pushing_self_cw)',
  'sd(hu_pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'ma[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)
# For indistinguishable Dyads
model_rownames_fixed_hu <- c(
    "Intercept", 
    "Hurdle Intercept",
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime",
    
    # HURDLE
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Hu Daily persuasion experienced", 
    "Hu Daily persuasion utilized (partner's view)", # OR partner received
    "Hu Daily pressure experienced", 
    "Hu Daily pressure utilized (partner's view)", 
    "Hu Daily pushing experienced", 
    "Hu Daily pushing utilized (partner's view)", 
    "Hu Day", 
    "Hu Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Hu Mean persuasion experienced", 
    "Hu Mean persuasion utilized (partner's view)", 
    "Hu Mean pressure experienced", 
    "Hu Mean pressure utilized (partner's view)", 
    "Hu Mean pushing experienced", 
    "Hu Mean pushing utilized (partner's view)", 
    "Hu Mean weartime"
  )


model_rownames_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Hurdle Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  
  # Hurdle
  "sd(Hu Daily persuasion experienced)", 
  "sd(Hu Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Hu Daily pressure experienced)", 
  "sd(Hu Daily pressure utilized (partner's view))", 
  "sd(Hu Daily pushing experienced)", 
  "sd(Hu Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'ma[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)
rows_to_pack_hu <- list(
  "Conditional Within-Person Effects" = c(3,10),
  "Conditional Between-Person Effects" = c(11,17),
  
  "Hurdle Within-Person Effects" = c(18,25),
  "Hurdle Between-Person Effects" = c(26,32),
  
  "Random Effects" = c(33, 46), 
  "Additional Parameters" = c(47,52)
  )

Subjective MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

Modelling using the gaussian family fails. Due to the many zeros, transformations won’t help estimating the models. We employ the negative binomial family.

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID) 
  , autocor = ~ arma(gr = coupleID, cov = TRUE)
  #, autocor = ~ unstr(gr = coupleID)
  #, autocor = ~ cosy(gr = coupleID)
)

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(4, 2)", class = "Intercept"), # for non-zero PA
  brms::set_prior("normal(6, 2.5)", class = "Intercept", dpar = 'hu'), # hurdle part
  
  brms::set_prior("normal(0, 5)", class = "sd", group = "coupleID", lb = 0), # Random effects SD for couple
  
  brms::set_prior("normal(0, 0.3)", class = "ar", lb = -1, ub = 1), # Autocorrelation prior
  brms::set_prior("normal(0, 0.3)", class = "ma", lb = -1, ub = 1), # moving average prior
  
  brms::set_prior("normal(0.1,2)", class = "sderr", lb = 0),
  brms::set_prior("normal(0.1,2)", class = "sigma", lb = 0) # Residual SD for lognormal
)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(),
  control = list(adapt_delta = 0.95),
  iter = 8000,
  warmup = 3000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "pa_sub_hu")
)
plot(pa_sub, ask = FALSE, nvariables = 2)

plot(pp_check(pa_sub, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(pa_sub))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

pp_check_transformed(pa_sub, transform = log1p)

loo_pa_sub <- loo(pa_sub)
loo_pa_sub
## 
## Computed from 20000 by 1342 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -5351.0  73.1
## p_loo       825.1  22.4
## looic     10702.1 146.1
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 1.6]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1058  78.8%   0       
##    (0.7, 1]   (bad)       271  20.2%   <NA>    
##    (1, Inf)   (very bad)   13   1.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  pa_sub, 
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T,
  include_p_direction = T) %>%
  print_df(rows_to_pack = rows_to_pack_hu
)
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 53.56* 43.52 65.98 1.000 9828.71 13100.94 1
Hurdle Intercept 0.34* 0.19 0.63 1.001 5899.52 8289.96 0.999
Conditional Within-Person Effects
Daily persuasion experienced 1.04* 1.00 1.08 1.000 13096.31 13973.69 0.983
Daily persuasion utilized (partner’s view) 1.04 1.00 1.08 1.000 13166.32 14722.86 0.966
Daily pressure experienced 0.94 0.84 1.04 1.000 12430.78 13130.58 0.885
Daily pressure utilized (partner’s view) 0.95 0.85 1.05 1.000 12806.25 13547.74 0.848
Daily pushing experienced 1.01 0.96 1.07 1.000 12288.43 13875.35 0.672
Daily pushing utilized (partner’s view) 1.00 0.95 1.06 1.000 12176.93 13067.96 0.565
Day 0.92 0.69 1.22 1.000 10805.99 13382.52 0.716
Daily weartime NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.08 0.73 1.57 1.000 8099.07 11263.07 0.654
Mean persuasion utilized (partner’s view) 1.11 0.76 1.61 1.000 7813.24 9463.86 0.706
Mean pressure experienced 0.88 0.40 1.91 1.000 8764.64 10070.47 0.631
Mean pressure utilized (partner’s view) 0.68 0.31 1.48 1.000 8586.42 10643.17 0.841
Mean pushing experienced 1.24 0.90 1.72 1.001 8116.63 10953.55 0.909
Mean pushing utilized (partner’s view) 1.22 0.88 1.68 1.000 8027.84 9704.30 0.886
Mean weartime NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 0.65* 0.50 0.81 1.000 18167.00 12325.72 1
Hu Daily persuasion utilized (partner’s view) 0.75* 0.59 0.92 1.000 15813.87 12981.63 0.997
Hu Daily pressure experienced 1.56 0.78 2.93 1.000 12350.77 9155.24 0.921
Hu Daily pressure utilized (partner’s view) 0.69 0.10 1.94 1.000 9690.05 8263.72 0.667
Hu Daily pushing experienced 0.67* 0.45 0.96 1.000 16658.28 13188.56 0.986
Hu Daily pushing utilized (partner’s view) 0.58* 0.36 0.84 1.000 14231.28 12560.84 0.998
Hu Day 1.27 0.75 2.17 1.001 31641.19 15366.82 0.814
Hu Daily weartime NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 0.64 0.13 3.12 1.001 4622.99 7947.19 0.716
Hu Mean persuasion utilized (partner’s view) 0.66 0.14 3.32 1.001 4643.77 7626.90 0.707
Hu Mean pressure experienced 4.19 0.62 37.41 1.000 7183.58 10013.58 0.924
Hu Mean pressure utilized (partner’s view) 3.45 0.52 31.05 1.000 7300.87 10348.27 0.895
Hu Mean pushing experienced 1.57 0.43 5.61 1.001 5817.77 9495.17 0.763
Hu Mean pushing utilized (partner’s view) 1.14 0.31 4.15 1.001 5694.21 9340.41 0.58
Hu Mean weartime NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.25 0.07 0.41 1.00 3596.94 3099.40 NA
sd(Hurdle Intercept) 1.47 1.06 2.04 1.00 6684.76 11223.26 NA
sd(Daily persuasion experienced) 0.03 0.00 0.08 1.00 4789.63 9208.46 NA
sd(Daily persuasion utilized (partner’s view)) 0.02 0.00 0.07 1.00 6239.09 9399.62 NA
sd(Daily pressure experienced) 0.06 0.00 0.20 1.00 8233.02 10928.12 NA
sd(Daily pressure utilized (partner’s view)) 0.07 0.00 0.21 1.00 7237.29 10383.03 NA
sd(Daily pushing experienced) 0.03 0.00 0.09 1.00 7047.09 9604.83 NA
sd(Daily pushing utilized (partner’s view)) 0.04 0.00 0.11 1.00 4575.13 9605.60 NA
sd(Hu Daily persuasion experienced) 0.26 0.01 0.61 1.00 5101.97 6854.84 NA
sd(Hu Daily persuasion utilized (partner’s view)) 0.25 0.01 0.55 1.00 6406.56 8018.65 NA
sd(Hu Daily pressure experienced) 0.53 0.02 1.80 1.00 7326.65 10118.74 NA
sd(Hu Daily pressure utilized (partner’s view)) 1.35 0.09 3.78 1.00 5814.81 5693.11 NA
sd(Hu Daily pushing experienced) 0.51 0.05 1.08 1.00 4919.13 6302.22 NA
sd(Hu Daily pushing utilized (partner’s view)) 0.55 0.05 1.24 1.00 4724.73 5976.66 NA
Additional Parameters
ar[1] 0.79 0.69 0.87 1.01 1913.99 5643.13 NA
ma[1] -0.09 -0.41 0.46 1.16 18.90 110.03 NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr 0.44 0.24 0.60 1.24 14.02 74.28 NA
sigma 0.33 0.10 0.47 1.25 13.60 54.74 NA

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

We tried negative binomial here as well for consistency, but the model did not converge. Poisson also did not work. As we have no zeros in this distribution, we log transform.

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = ~ arma(gr = coupleID, cov = TRUE)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(4, 2)", class = "Intercept"), # for non-zero PA
  #brms::set_prior("normal(0, 2.5)", class = "Intercept", dpar = 'hu'), # hurdle part
  
  brms::set_prior("normal(0, 5)", class = "sd", group = "coupleID", lb = 0), # Random effects SD for couple
  
  brms::set_prior("normal(0, 0.3)", class = "ar", lb = -1, ub = 1), # Autocorrelation prior
  brms::set_prior("normal(0, 0.3)", class = "ma", lb = -1, ub = 1), # moving average prior
  
  brms::set_prior("normal(0.1,2)", class = "sderr", lb = 0),
  brms::set_prior("student_t(3, 0, 10)", class = "sigma", lb = 0) # Residual SD for lognormal
)


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::lognormal(),
  control = list(adapt_delta = 0.95),
  iter = 8000,
  warmup = 3000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "pa_obj_lognormal")
  #, file_refit = 'always'
)
plot(pa_obj_log, ask = FALSE, nvariables = 2)

plot(pp_check(pa_obj_log, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(pa_obj_log))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(pa_obj_log, transform = log1p)

loo_pa_obj_log <- loo(pa_obj_log)
loo_pa_obj_log
## 
## Computed from 12000 by 1214 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -6849.9 39.0
## p_loo       401.9 19.3
## looic     13699.8 78.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 0.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1190  98.0%   1       
##    (0.7, 1]   (bad)        23   1.9%   <NA>    
##    (1, Inf)   (very bad)    1   0.1%   <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  pa_obj_log, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 123.11* 107.69 140.76 1.000 4965.60 6025.57
Within-Person Effects
Daily persuasion experienced 1.01 0.98 1.05 1.000 8215.74 7621.04
Daily persuasion utilized (partner’s view) 1.02 0.98 1.06 1.000 10555.20 8230.76
Daily pressure experienced 0.97 0.87 1.07 1.001 9625.90 7468.88
Daily pressure utilized (partner’s view) 0.94 0.86 1.03 1.000 9492.95 8064.57
Daily pushing experienced 1.03 0.97 1.09 1.000 9025.58 7734.85
Daily pushing utilized (partner’s view) 1.02 0.97 1.08 1.000 10386.80 9285.70
Day 0.97 0.83 1.13 1.000 10311.21 9047.61
Daily weartime 1.00* 1.00 1.00 1.001 10175.12 9143.71
Between-Person Effects
Mean persuasion experienced 1.05 0.77 1.45 1.000 3600.09 5551.97
Mean persuasion utilized (partner’s view) 0.87 0.63 1.19 1.000 3538.01 5845.54
Mean pressure experienced 0.78 0.55 1.13 1.001 5743.15 8068.77
Mean pressure utilized (partner’s view) 1.01 0.72 1.45 1.001 5185.46 7697.84
Mean pushing experienced 1.10 0.84 1.44 1.002 4042.39 6190.50
Mean pushing utilized (partner’s view) 1.18 0.91 1.55 1.002 3890.84 5577.78
Mean weartime 1.00 1.00 1.00 1.000 12630.10 8837.64
Random Effects
sd(Intercept) 0.29 0.18 0.40 1.00 1982.96 1920.07
sd(Daily persuasion experienced) 0.05 0.00 0.10 1.00 2776.36 3295.17
sd(Daily persuasion utilized (partner’s view)) 0.04 0.00 0.10 1.00 3009.58 4583.66
sd(Daily pressure experienced) 0.08 0.00 0.23 1.00 4202.78 6703.78
sd(Daily pressure utilized (partner’s view)) 0.05 0.00 0.14 1.00 6680.95 5730.26
sd(Daily pushing experienced) 0.06 0.00 0.14 1.00 2864.78 5405.82
sd(Daily pushing utilized (partner’s view)) 0.04 0.00 0.10 1.00 4521.24 5729.57
Additional Parameters
ar[1] 0.75 0.54 0.91 1.02 646.14 1354.20
ma[1] -0.12 -0.60 0.51 1.06 54.35 104.05
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr 0.24 0.10 0.47 1.10 31.29 30.36
sigma 0.44 0.25 0.51 1.09 32.69 30.75

Affect

range(df_double$aff, na.rm = T) 
## [1] 1 6
hist(df_double$aff, breaks = 15)

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = ~ arma(gr = coupleID, cov = TRUE)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("normal(0, 0.3)", class = "ar", lb = -1, ub = 1), # Autocorrelation prior
  brms::set_prior("normal(0, 0.3)", class = "ma", lb = -1, ub = 1), # moving average prior
  
  brms::set_prior("normal(0.1,2)", class = "sderr", lb = 0),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::skew_normal(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "mood_skew")
  #, file_refit = 'always'
)
plot(mood, ask = FALSE, nvariables = 2)

plot(pp_check(mood, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(mood))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(mood, transform = log1p)

loo_mood <- loo(mood)
loo_mood
## 
## Computed from 12000 by 1342 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1444.9 35.8
## p_loo        65.7  6.7
## looic      2889.8 71.6
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.3]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1264  94.2%   706     
##    (0.7, 1]   (bad)        71   5.3%   <NA>    
##    (1, Inf)   (very bad)    7   0.5%   <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  mood, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.89* 4.74 5.03 1.004 1573.16 2579.54
Within-Person Effects
Daily persuasion experienced 0.00 -0.03 0.03 1.000 11049.26 8352.67
Daily persuasion utilized (partner’s view) 0.00 -0.02 0.03 1.000 11245.39 8678.10
Daily pressure experienced -0.01 -0.11 0.10 1.000 7427.04 8066.53
Daily pressure utilized (partner’s view) -0.09 -0.23 0.03 1.001 5648.83 7576.12
Daily pushing experienced 0.02 -0.02 0.07 1.001 9623.09 8776.16
Daily pushing utilized (partner’s view) 0.05* 0.00 0.10 1.000 7798.79 7727.69
Day -0.01 -0.10 0.09 1.000 12064.82 8720.80
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.23 -0.22 0.67 1.001 2439.73 4065.52
Mean persuasion utilized (partner’s view) -0.02 -0.47 0.43 1.001 2335.84 3881.05
Mean pressure experienced -0.15 -0.72 0.41 1.002 2993.16 4716.78
Mean pressure utilized (partner’s view) 0.03 -0.53 0.59 1.001 3235.87 4645.00
Mean pushing experienced -0.02 -0.38 0.35 1.004 2787.59 4569.35
Mean pushing utilized (partner’s view) 0.08 -0.29 0.45 1.002 2620.47 4828.78
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.41 0.31 0.53 1.00 3732.68 6036.53
sd(Daily persuasion experienced) 0.02 0.00 0.05 1.00 4953.25 5663.82
sd(Daily persuasion utilized (partner’s view)) 0.02 0.00 0.04 1.00 6483.28 6724.47
sd(Daily pressure experienced) 0.08 0.00 0.24 1.00 5003.45 6791.59
sd(Daily pressure utilized (partner’s view)) 0.12 0.01 0.30 1.00 4296.09 4747.85
sd(Daily pushing experienced) 0.04 0.00 0.10 1.00 4254.31 4964.34
sd(Daily pushing utilized (partner’s view)) 0.07 0.02 0.13 1.00 6584.47 5671.86
Additional Parameters
ar[1] -1.00 -1.00 -0.99 1.00 3133.25 4808.19
ma[1] 0.04 -0.52 0.54 1.00 14102.68 7613.32
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr 0.01 0.00 0.03 1.00 3803.60 5904.53
sigma 0.80 0.76 0.83 1.00 12164.25 8750.51

Reactance

Gaussian

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 6)

formula <- bf(
  reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID), 
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = ~ arma(gr = coupleID, cov = TRUE)
)



prior1 <- c(
  brms::set_prior("normal(0, 1.5)", class = "b"),
  brms::set_prior("normal(0, 2)", class = "Intercept"), # for non-zero PA
  brms::set_prior("normal(40, 10)", class = "Intercept", dpar = 'hu'), # hurdle part
  
  brms::set_prior("normal(0, 5)", class = "sd", group = "coupleID", lb = 0), # Random effects SD for couple
  
  brms::set_prior("normal(0, 0.3)", class = "ar", lb = -1, ub = 1), # Autocorrelation prior
  brms::set_prior("normal(0, 0.3)", class = "ma", lb = -1, ub = 1), # moving average prior
  
  brms::set_prior("normal(0.1,2)", class = "sderr", lb = 0),
  brms::set_prior("normal(0,5)", class = "sigma", lb = 0) # Residual SD for lognormal
)




#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "reactance_hu")
  #, file_refit = 'always'
)
plot(reactance, ask = FALSE, nvariables = 2)

plot(pp_check(reactance, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(reactance))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

pp_check_transformed(reactance, transform = log1p)

loo_reactance <- loo(reactance)
loo_reactance
## 
## Computed from 12000 by 403 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -350.1 23.8
## p_loo       161.7 12.6
## looic       700.3 47.6
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 0.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     335   83.1%   1       
##    (0.7, 1]   (bad)       65   16.1%   <NA>    
##    (1, Inf)   (very bad)   3    0.7%   <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  reactance, 
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_hu)
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.96* 1.50 2.54 1.001 3892.12 8347.45
Hurdle Intercept 5.42* 1.94 16.61 1.004 2012.22 5131.61
Conditional Within-Person Effects
Daily persuasion experienced 0.99 0.89 1.12 1.001 3158.00 6658.63
Daily persuasion utilized (partner’s view) 1.00 0.88 1.13 1.002 1431.21 705.28
Daily pressure experienced 1.00 0.84 1.18 1.009 582.80 269.42
Daily pressure utilized (partner’s view) 0.99 0.64 1.54 1.001 2656.72 4092.91
Daily pushing experienced 1.03 0.92 1.15 1.002 3812.68 6919.63
Daily pushing utilized (partner’s view) 1.01 0.82 1.25 1.001 4085.63 7374.25
Day 0.99 0.64 1.54 1.001 5041.26 8074.46
Daily weartime NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 0.74 0.43 1.26 1.003 3149.35 6142.32
Mean persuasion utilized (partner’s view) 1.09 0.66 1.83 1.001 4049.83 7279.90
Mean pressure experienced 1.47 0.88 2.47 1.004 3540.08 6764.35
Mean pressure utilized (partner’s view) 1.09 0.65 1.81 1.003 1783.32 3710.65
Mean pushing experienced 0.89 0.65 1.23 1.001 2561.54 8141.97
Mean pushing utilized (partner’s view) 1.01 0.72 1.41 1.002 3782.31 7609.55
Mean weartime NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 1.56* 1.11 2.32 1.002 1946.52 3701.37
Hu Daily persuasion utilized (partner’s view) 1.04 0.66 1.79 1.001 5306.92 6772.91
Hu Daily pressure experienced 0.57 0.18 1.76 1.003 1622.79 2968.22
Hu Daily pressure utilized (partner’s view) 1.00 0.12 13.09 1.015 326.58 99.07
Hu Daily pushing experienced 0.72 0.49 1.04 1.004 870.71 2246.16
Hu Daily pushing utilized (partner’s view) 1.20 0.72 2.03 1.006 934.68 2640.79
Hu Day 0.75 0.22 2.53 1.002 3550.56 6253.59
Hu Daily weartime NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 0.63 0.07 4.98 1.002 2576.11 6427.43
Hu Mean persuasion utilized (partner’s view) 0.72 0.07 6.88 1.005 903.84 285.74
Hu Mean pressure experienced 0.05 0.00 1.28 1.001 2500.11 3294.38
Hu Mean pressure utilized (partner’s view) 3.28 0.14 157.82 1.003 2255.32 6120.90
Hu Mean pushing experienced 0.80 0.16 5.31 1.006 650.29 126.44
Hu Mean pushing utilized (partner’s view) 1.28 0.16 8.89 1.008 460.85 264.71
Hu Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.13 0.01 0.35 1.00 2965.17 5684.96
sd(Hurdle Intercept) 1.88 0.96 3.48 1.01 419.61 104.48
sd(Daily persuasion experienced) 0.15 0.02 0.29 1.01 1742.95 2671.81
sd(Daily persuasion utilized (partner’s view)) 0.06 0.00 0.19 1.00 1501.92 4975.99
sd(Daily pressure experienced) 0.10 0.00 0.31 1.00 3998.43 5767.26
sd(Daily pressure utilized (partner’s view)) 0.28 0.01 1.11 1.00 3686.85 3436.33
sd(Daily pushing experienced) 0.09 0.00 0.24 1.00 2475.90 5946.62
sd(Daily pushing utilized (partner’s view)) 0.22 0.03 0.45 1.01 1123.99 3348.94
sd(Hu Daily persuasion experienced) 0.45 0.05 0.93 1.00 2503.92 3876.32
sd(Hu Daily persuasion utilized (partner’s view)) 0.70 0.08 1.45 1.00 1987.62 1455.29
sd(Hu Daily pressure experienced) 1.50 0.25 3.51 1.00 2856.04 3863.46
sd(Hu Daily pressure utilized (partner’s view)) 1.59 0.05 4.98 1.01 615.86 264.17
sd(Hu Daily pushing experienced) 0.43 0.03 1.00 1.00 1159.22 3617.78
sd(Hu Daily pushing utilized (partner’s view)) 0.41 0.02 1.17 1.00 929.46 259.25
Additional Parameters
ar[1] 0.08 -0.42 0.54 1.00 1580.35 1389.03
ma[1] 0.10 -0.46 0.62 1.00 2969.54 4465.57
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr 0.25 0.01 0.45 1.02 96.30 169.93
sigma 0.31 0.06 0.49 1.05 63.97 17.31 
hypothesis(reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... > 0    -0.03      0.12    -0.23     0.17       0.72
##   Post.Prob Star
## 1      0.42     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(reactance, "hu_pressure_self_cw > hu_pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (hu_pressure_self... > 0    -0.23      0.62    -1.22     0.74        0.5
##   Post.Prob Star
## 1      0.33     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = ~ arma(gr = coupleID, cov = TRUE)
)



prior1 <- c(
  set_prior("normal(0, 2.5)", class = "b"),
  set_prior("normal(0, 5)", class = "Intercept"),
  set_prior("normal(0, 1)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("normal(0, 0.3)", class = "ar", lb = -1, ub = 1), # Autocorrelation prior
  brms::set_prior("normal(0, 0.3)", class = "ma", lb = -1, ub = 1), # moving average prior
  
  brms::set_prior("normal(0.1,2)", class = "sderr", lb = 0)
)




#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  control = list(adapt_delta = 0.95),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "reactance_ordinal_noar")
  #, file_refit = 'always'
)
plot(reactance_ordinal, ask = FALSE, nvariables = 2)

plot(pp_check(reactance_ordinal, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(reactance_ordinal))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(reactance_ordinal, transform = log1p)

loo_reactance_ordinal <- loo(reactance_ordinal)
loo_reactance_ordinal
## 
## Computed from 12000 by 403 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -344.9 22.9
## p_loo       180.4 14.0
## looic       689.8 45.7
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     198   49.1%   104     
##    (0.7, 1]   (bad)      199   49.4%   <NA>    
##    (1, Inf)   (very bad)   6    1.5%   <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  reactance_ordinal, 
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercepts
Intercept NA NA NA NA NA NA
Intercept[1] 12.57* 2.72 80.47 1.000 2607.69 5015.89
Intercept[2] 49.58* 8.15 458.69 1.001 1578.09 3842.69
Intercept[3] 300.44* 30.86 4802.89 1.001 1180.99 2068.45
Intercept[4] 3418.94* 174.50 121990.99 1.001 1016.54 1736.99
Intercept[5] 140267.36* 2501.00 16508154.02 1.001 1048.96 1941.23
Within-Person Effects
Daily persuasion experienced 0.48* 0.27 0.77 1.001 3245.33 6087.86
Daily persuasion utilized (partner’s view) 1.01 0.54 1.81 1.001 6705.40 6276.08
Daily pressure experienced 2.62 0.94 7.44 1.000 3551.40 6632.25
Daily pressure utilized (partner’s view) 1.53 0.41 5.23 1.000 7930.81 8699.38
Daily pushing experienced 1.62 0.91 3.09 1.000 4741.04 6378.06
Daily pushing utilized (partner’s view) 0.75 0.36 1.51 1.000 8188.03 7066.32
Day 1.29 0.15 11.32 1.000 9907.70 9228.37
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.01 0.11 10.29 1.001 6994.12 7429.84
Mean persuasion utilized (partner’s view) 1.24 0.11 14.84 1.000 7513.11 8129.01
Mean pressure experienced 3.12 0.26 43.70 1.001 7507.68 8393.24
Mean pressure utilized (partner’s view) 1.24 0.09 15.81 1.000 8110.18 8245.48
Mean pushing experienced 1.91 0.33 12.73 1.001 6837.67 6440.44
Mean pushing utilized (partner’s view) 0.87 0.12 6.69 1.001 7303.88 7532.68
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.41 0.25 2.52 1.00 2215.09 2571.70
sd(Daily persuasion experienced) 0.34 0.01 0.93 1.00 2792.32 4880.69
sd(Daily persuasion utilized (partner’s view)) 0.48 0.02 1.24 1.00 2772.51 5526.35
sd(Daily pressure experienced) 0.90 0.06 2.07 1.00 2908.71 5159.82
sd(Daily pressure utilized (partner’s view)) 0.72 0.03 1.95 1.00 5764.57 5991.92
sd(Daily pushing experienced) 0.63 0.05 1.42 1.00 2476.84 4454.31
sd(Daily pushing utilized (partner’s view)) 0.34 0.01 1.01 1.00 6069.02 6393.85
Additional Parameters
ar[1] 0.28 -0.21 0.71 1.00 2648.16 3971.89
ma[1] 0.05 -0.36 0.52 1.00 4385.44 7031.10
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr 3.23 1.37 5.09 1.00 764.01 1243.91
sigma NA NA NA NA NA NA
disc 1.00 1.00 1.00 NA NA NA 
hypothesis(reactance_ordinal, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... > 0     0.48      0.63    -0.56     1.53       3.63
##   Post.Prob Star
## 1      0.78     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = ~ ar(gr = coupleID)
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1)
)


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "is_reactance_noar")
  #, file_refit = 'always'
)
plot(is_reactance, ask = FALSE, nvariables = 2)

plot(pp_check(is_reactance, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(is_reactance))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(is_reactance, transform = log1p)

loo_is_reactance <- loo(is_reactance)
loo_is_reactance
## 
## Computed from 12000 by 403 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -367.2 19.2
## p_loo       323.6 17.5
## looic       734.3 38.3
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)       5    1.2%   749     
##    (0.7, 1]   (bad)      145   36.0%   <NA>    
##    (1, Inf)   (very bad) 253   62.8%   <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  is_reactance, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
OR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.33 0.04 2.42 1.000 8136.75 7710.30
Within-Person Effects
Daily persuasion experienced 0.30* 0.10 0.72 1.000 5379.22 7183.33
Daily persuasion utilized (partner’s view) 1.72 0.53 7.52 1.001 5921.25 6811.67
Daily pressure experienced 5.74 0.66 65.46 1.001 6505.22 7514.46
Daily pressure utilized (partner’s view) 1.61 0.15 18.88 1.001 8701.15 8623.44
Daily pushing experienced 2.96* 1.08 11.61 1.001 5036.37 6193.62
Daily pushing utilized (partner’s view) 0.77 0.20 3.52 1.000 6747.95 6865.75
Day 1.99 0.10 40.64 1.000 9224.50 8920.00
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 8.01 0.29 226.79 1.000 8074.54 8987.13
Mean persuasion utilized (partner’s view) 4.10 0.14 122.71 1.001 9690.86 8710.65
Mean pressure experienced 15.71 0.27 989.74 1.000 11828.56 9544.85
Mean pressure utilized (partner’s view) 0.75 0.01 53.46 1.000 14614.32 9626.17
Mean pushing experienced 10.20 0.58 209.93 1.000 7280.14 8477.24
Mean pushing utilized (partner’s view) 1.86 0.08 41.76 1.000 8290.11 8913.90
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 4.76 2.89 6.98 1.00 4027.89 6109.26
sd(Daily persuasion experienced) 0.81 0.04 1.99 1.00 2681.73 4564.17
sd(Daily persuasion utilized (partner’s view)) 1.53 0.18 3.36 1.00 3233.18 3914.65
sd(Daily pressure experienced) 2.40 0.26 4.94 1.00 3293.43 2816.17
sd(Daily pressure utilized (partner’s view)) 1.61 0.06 4.33 1.00 5179.41 5307.73
sd(Daily pushing experienced) 1.01 0.06 2.49 1.00 2565.30 4926.10
sd(Daily pushing utilized (partner’s view)) 0.91 0.04 2.57 1.00 4262.83 6155.71
Additional Parameters
ar[1] 0.26 -0.07 0.78 1.00 1547.00 1167.99
ma[1] NA NA NA NA NA NA
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr 5.42 2.74 8.78 1.00 2280.05 3184.86
sigma NA NA NA NA NA NA 
hypothesis(is_reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... > 0     0.66      1.34    -1.51     2.88       2.36
##   Post.Prob Star
## 1       0.7     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

if (report_ordinal) {
  model_rows_random_final <- model_rows_random_ordinal
  model_rows_fixed_final <- model_rows_fixed_ordinal
  model_rownames_fixed_final <- model_rownames_fixed_ordinal
  model_rownames_random_final <- model_rownames_random_ordinal
  rows_to_pack_final <- rows_to_pack_ordinal
} else {
  model_rows_random_final <- model_rows_random_hu
  model_rows_fixed_final <- model_rows_fixed_hu
  model_rownames_fixed_final <- model_rownames_fixed_hu
  model_rownames_random_final <- model_rownames_random_hu
  rows_to_pack_final <- rows_to_pack_hu
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood,
  reactance,
  reactance_ordinal,
  is_reactance,
  
  model_rows_random = model_rows_random_final,
  model_rows_fixed = model_rows_fixed_final,
  model_rownames_random = model_rownames_random_final,
  model_rownames_fixed = model_rownames_fixed_final
) 
## [1] "pa_sub"
## [1] "pa_obj_log"
## [1] "mood"
## [1] "reactance"
## [1] "reactance_ordinal"
## [1] "is_reactance"
# pretty printing
summary_all_models <- all_models %>%
  print_df(rows_to_pack = rows_to_pack_final) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = 2, 
      "Device-Based MVPA Lognormal" = 2, 
      "Mood Skewnormal" = 2,
      "Reactance Lognormal" = 2, 
      "Reactance Ordinal" = 2,
      "Reactance Dichotome" = 2)
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack_final,
  file.path("Output", "AllModels_experimental.xlsx"), 
  merge_option = 'both', 
  simplify_2nd_row = TRUE,
  colwidths = c(38, 7.2, 13.3, 7.2, 13.3, 7.2, 13.3,7.2, 13.3,7.2, 13.3,7.2, 13.3),
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
## 
##     guess_encoding
summary_all_models
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Lognormal
Mood Skewnormal
Reactance Lognormal
Reactance Ordinal
Reactance Dichotome
exp(Est.) pa_sub 95% CI pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log b mood 95% CI mood exp(Est.) reactance 95% CI reactance OR reactance_ordinal 95% CI reactance_ordinal OR is_reactance 95% CI is_reactance
Intercept 53.56* [43.52, 65.98] 123.11* [107.69, 140.76] 4.89* [ 4.74, 5.03] 1.96* [1.50, 2.54] NA NA 0.33 [0.04, 2.42]
Hurdle Intercept 0.34* [ 0.19, 0.63] NA NA NA NA 5.42* [1.94, 16.61] NA NA NA NA
Conditional Within-Person Effects
Daily persuasion experienced 1.04* [ 1.00, 1.08] 1.01 [ 0.98, 1.05] 0.00 [-0.03, 0.03] 0.99 [0.89, 1.12] 0.48* [0.27, 0.77] 0.30* [0.10, 0.72]
Daily persuasion utilized (partner’s view) 1.04 [ 1.00, 1.08] 1.02 [ 0.98, 1.06] 0.00 [-0.02, 0.03] 1.00 [0.88, 1.13] 1.01 [0.54, 1.81] 1.72 [0.53, 7.52]
Daily pressure experienced 0.94 [ 0.84, 1.04] 0.97 [ 0.87, 1.07] -0.01 [-0.11, 0.10] 1.00 [0.84, 1.18] 2.62 [0.94, 7.44] 5.74 [0.66, 65.46]
Daily pressure utilized (partner’s view) 0.95 [ 0.85, 1.05] 0.94 [ 0.86, 1.03] -0.09 [-0.23, 0.03] 0.99 [0.64, 1.54] 1.53 [0.41, 5.23] 1.61 [0.15, 18.88]
Daily pushing experienced 1.01 [ 0.96, 1.07] 1.03 [ 0.97, 1.09] 0.02 [-0.02, 0.07] 1.03 [0.92, 1.15] 1.62 [0.91, 3.09] 2.96* [1.08, 11.61]
Daily pushing utilized (partner’s view) 1.00 [ 0.95, 1.06] 1.02 [ 0.97, 1.08] 0.05* [ 0.00, 0.10] 1.01 [0.82, 1.25] 0.75 [0.36, 1.51] 0.77 [0.20, 3.52]
Day 0.92 [ 0.69, 1.22] 0.97 [ 0.83, 1.13] -0.01 [-0.10, 0.09] 0.99 [0.64, 1.54] 1.29 [0.15, 11.32] 1.99 [0.10, 40.64]
Daily weartime NA NA 1.00* [ 1.00, 1.00] NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.08 [ 0.73, 1.57] 1.05 [ 0.77, 1.45] 0.23 [-0.22, 0.67] 0.74 [0.43, 1.26] 1.01 [0.11, 10.29] 8.01 [0.29, 226.79]
Mean persuasion utilized (partner’s view) 1.11 [ 0.76, 1.61] 0.87 [ 0.63, 1.19] -0.02 [-0.47, 0.43] 1.09 [0.66, 1.83] 1.24 [0.11, 14.84] 4.10 [0.14, 122.71]
Mean pressure experienced 0.88 [ 0.40, 1.91] 0.78 [ 0.55, 1.13] -0.15 [-0.72, 0.41] 1.47 [0.88, 2.47] 3.12 [0.26, 43.70] 15.71 [0.27, 989.74]
Mean pressure utilized (partner’s view) 0.68 [ 0.31, 1.48] 1.01 [ 0.72, 1.45] 0.03 [-0.53, 0.59] 1.09 [0.65, 1.81] 1.24 [0.09, 15.81] 0.75 [0.01, 53.46]
Mean pushing experienced 1.24 [ 0.90, 1.72] 1.10 [ 0.84, 1.44] -0.02 [-0.38, 0.35] 0.89 [0.65, 1.23] 1.91 [0.33, 12.73] 10.20 [0.58, 209.93]
Mean pushing utilized (partner’s view) 1.22 [ 0.88, 1.68] 1.18 [ 0.91, 1.55] 0.08 [-0.29, 0.45] 1.01 [0.72, 1.41] 0.87 [0.12, 6.69] 1.86 [0.08, 41.76]
Mean weartime NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 0.65* [ 0.50, 0.81] NA NA NA NA 1.56* [1.11, 2.32] NA NA NA NA
Hu Daily persuasion utilized (partner’s view) 0.75* [ 0.59, 0.92] NA NA NA NA 1.04 [0.66, 1.79] NA NA NA NA
Hu Daily pressure experienced 1.56 [ 0.78, 2.93] NA NA NA NA 0.57 [0.18, 1.76] NA NA NA NA
Hu Daily pressure utilized (partner’s view) 0.69 [ 0.10, 1.94] NA NA NA NA 1.00 [0.12, 13.09] NA NA NA NA
Hu Daily pushing experienced 0.67* [ 0.45, 0.96] NA NA NA NA 0.72 [0.49, 1.04] NA NA NA NA
Hu Daily pushing utilized (partner’s view) 0.58* [ 0.36, 0.84] NA NA NA NA 1.20 [0.72, 2.03] NA NA NA NA
Hu Day 1.27 [ 0.75, 2.17] NA NA NA NA 0.75 [0.22, 2.53] NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 0.64 [ 0.13, 3.12] NA NA NA NA 0.63 [0.07, 4.98] NA NA NA NA
Hu Mean persuasion utilized (partner’s view) 0.66 [ 0.14, 3.32] NA NA NA NA 0.72 [0.07, 6.88] NA NA NA NA
Hu Mean pressure experienced 4.19 [ 0.62, 37.41] NA NA NA NA 0.05 [0.00, 1.28] NA NA NA NA
Hu Mean pressure utilized (partner’s view) 3.45 [ 0.52, 31.05] NA NA NA NA 3.28 [0.14, 157.82] NA NA NA NA
Hu Mean pushing experienced 1.57 [ 0.43, 5.61] NA NA NA NA 0.80 [0.16, 5.31] NA NA NA NA
Hu Mean pushing utilized (partner’s view) 1.14 [ 0.31, 4.15] NA NA NA NA 1.28 [0.16, 8.89] NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.25 [ 0.07, 0.41] 0.29 [ 0.18, 0.40] 0.41 [ 0.31, 0.53] 0.13 [ 0.01, 0.35] 1.41 [ 0.25, 2.52] 4.76 [ 2.89, 6.98]
sd(Hurdle Intercept) 1.47 [ 1.06, 2.04] NA NA NA NA 1.88 [ 0.96, 3.48] NA NA NA NA
sd(Daily persuasion experienced) 0.03 [ 0.00, 0.08] 0.05 [ 0.00, 0.10] 0.02 [ 0.00, 0.05] 0.15 [ 0.02, 0.29] 0.34 [ 0.01, 0.93] 0.81 [ 0.04, 1.99]
sd(Daily persuasion utilized (partner’s view)) 0.02 [ 0.00, 0.07] 0.04 [ 0.00, 0.10] 0.02 [ 0.00, 0.04] 0.06 [ 0.00, 0.19] 0.48 [ 0.02, 1.24] 1.53 [ 0.18, 3.36]
sd(Daily pressure experienced) 0.06 [ 0.00, 0.20] 0.08 [ 0.00, 0.23] 0.08 [ 0.00, 0.24] 0.10 [ 0.00, 0.31] 0.90 [ 0.06, 2.07] 2.40 [ 0.26, 4.94]
sd(Daily pressure utilized (partner’s view)) 0.07 [ 0.00, 0.21] 0.05 [ 0.00, 0.14] 0.12 [ 0.01, 0.30] 0.28 [ 0.01, 1.11] 0.72 [ 0.03, 1.95] 1.61 [ 0.06, 4.33]
sd(Daily pushing experienced) 0.03 [ 0.00, 0.09] 0.06 [ 0.00, 0.14] 0.04 [ 0.00, 0.10] 0.09 [ 0.00, 0.24] 0.63 [ 0.05, 1.42] 1.01 [ 0.06, 2.49]
sd(Daily pushing utilized (partner’s view)) 0.04 [ 0.00, 0.11] 0.04 [ 0.00, 0.10] 0.07 [ 0.02, 0.13] 0.22 [ 0.03, 0.45] 0.34 [ 0.01, 1.01] 0.91 [ 0.04, 2.57]
sd(Hu Daily persuasion experienced) 0.26 [ 0.01, 0.61] NA NA NA NA 0.45 [ 0.05, 0.93] NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) 0.25 [ 0.01, 0.55] NA NA NA NA 0.70 [ 0.08, 1.45] NA NA NA NA
sd(Hu Daily pressure experienced) 0.53 [ 0.02, 1.80] NA NA NA NA 1.50 [ 0.25, 3.51] NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) 1.35 [ 0.09, 3.78] NA NA NA NA 1.59 [ 0.05, 4.98] NA NA NA NA
sd(Hu Daily pushing experienced) 0.51 [ 0.05, 1.08] NA NA NA NA 0.43 [ 0.03, 1.00] NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) 0.55 [ 0.05, 1.24] NA NA NA NA 0.41 [ 0.02, 1.17] NA NA NA NA
Additional Parameters
ar[1] 0.79 [ 0.69, 0.87] 0.75 [ 0.54, 0.91] -1.00 [-1.00, -0.99] 0.08 [-0.42, 0.54] 0.28 [-0.21, 0.71] 0.26 [-0.07, 0.78]
ma[1] -0.09 [-0.41, 0.46] -0.12 [-0.60, 0.51] 0.04 [-0.52, 0.54] 0.10 [-0.46, 0.62] 0.05 [-0.36, 0.52] NA NA
nu NA NA NA NA NA NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA NA NA NA NA NA
sderr 0.44 [ 0.24, 0.60] 0.24 [ 0.10, 0.47] 0.01 [ 0.00, 0.03] 0.25 [ 0.01, 0.45] 3.23 [ 1.37, 5.09] 5.42 [ 2.74, 8.78]
sigma 0.33 [ 0.10, 0.47] 0.44 [ 0.25, 0.51] 0.80 [ 0.76, 0.83] 0.31 [ 0.06, 0.49] NA NA NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::cite_packages()